Results for Hitchcock & Frank 2023
Peter Frank Hitchcock - LNCC Lab, Brown University
2024-01-24
- Load data from studies 1 and 2
- Behavioral Results for Learning and Test phase
- Modeling
Results (Note: model numbers match those from
s.Rbut not those used in paper…)- Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter
- Model 3 — Same as m1 but not letting decay vary
- Model 4 — Same as m3 but to 0 inits
- Model 5 — Same as m1 but lets beta vary instead of decay
- Model 6 — Same as m1 but lets learning rate vary instead of decay
- Model 11 — Simple q-learner model
- Model 12 — Simple Bayes model
- Model 13 — M12 but allowing beta to vary
- Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond
- Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond
- Model 27 — Fixing both ES and epsilon at the group level
- Model 29 — Model 1 but with no choice kernel
- Model Comparison
sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"doParallel",
"patchwork",
"lme4",
"lmerTest",
"testit",
"ggpubr",
"latex2exp"
),
require, character=TRUE)## Loading required package: rjson## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: ggplot2## Loading required package: stringr## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: doParallel## Loading required package: iterators## Loading required package: parallel## Loading required package: patchwork## Loading required package: lme4## Loading required package: Matrix## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: testit## Loading required package: ggpubr## Loading required package: latex2exp## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## doParallel patchwork lme4 lmerTest testit ggpubr latex2exp
## TRUE TRUE TRUE TRUE TRUE TRUE TRUEsf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()## ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value ? ?
## visible FALSE FALSE
## ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value ? ?
## visible FALSE FALSE
## ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearner.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/aModelUtils.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModel.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/RecodeDfs.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayQVals.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value ?
## visible FALSEDefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))Load data from studies 1 and 2
s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
#read.csv("../data/cleaned_files/s1_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv") #read.csv("../data/cleaned_files/s2_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv") Harmonize variables and create some separate vars
# Study 2 harmonize
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive"
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2
s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))
s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))Define paths and functions
# MLE results
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations (old folder name — includes emp bayes res)
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function
rm <- function(path, model_str) read.csv(paste0(path, model_str))Behavioral Results for Learning and Test phase
Visualization
Summarize training mean and within-subject SEM
s1_train_summs <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_tr_summs <- Rmisc::summarySEwithin(s1_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_train_summs <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions1_test_summs <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_te_summs <- Rmisc::summarySEwithin(s1_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_test_summs <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions1_tr_summs$experiment <- 1
s1_tr_summs$phase <- "Learning"
s2_tr_summs$experiment <- 2
s2_tr_summs$phase <- "Learning"
s1_te_summs$experiment <- 1
s1_te_summs$phase <- "Test"
s2_te_summs$experiment <- 2
s2_te_summs$phase <- "Test"
comb_summs <- rbind(s1_tr_summs, s2_tr_summs)
comb_summs$experiment <- factor(comb_summs$experiment)
comb_te_summs <- rbind(s1_te_summs, s2_te_summs)
comb_te_summs$experiment <- factor(comb_te_summs$experiment)a <- ggplot(comb_summs,
aes(x=experiment, y=m, group=condition, fill=condition)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) +
geom_hline(yintercept=.5, size=2) +
scale_fill_manual(values=c("purple", "orange")) +
geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
group=condition),
inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
geom_point(size=4, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol +
xlab("") + ylab("Proportion correct") + ggtitle("Learning") + tp +
ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.a
b <- ggplot(comb_te_summs,
aes(x=experiment, y=m, group=condition, fill=condition)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) +
geom_hline(yintercept=.5, size=2) +
scale_fill_manual(values=c("purple", "orange")) +
geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
group=condition),
inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
geom_point(size=3, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol + tp +
#lp + Used to cut legend and then turned off so size is comparable
xlab("") + ylab("") + ggtitle("Test") +
ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
b
perf_comb <- a + bperf_comb
#ggsave("../paper/figs/fig_parts/perf_plot.png", perf_comb, height=6, width=10, dpi=700)For resub, showing learning curves in Fig 2
s1_train_summs_lc <- s1_train %>% group_by(condition, ID, trial_within_condition, probability) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'ID',
## 'trial_within_condition'. You can override using the `.groups` argument.s1_tr_summs_lc <- Rmisc::summarySEwithin(s1_train_summs_lc,
measurevar = "m",
withinvars = c("condition", "trial_within_condition", "probability"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, trial_within_condition, probabilitys2_train_summs_lc <- s2_train %>% group_by(condition, ID, trial_within_condition, probability) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition', 'ID',
## 'trial_within_condition'. You can override using the `.groups` argument.s2_tr_summs_lc <- Rmisc::summarySEwithin(s2_train_summs_lc,
measurevar = "m",
withinvars = c("condition", "trial_within_condition", "probability"),
idvar = "ID")## Automatically converting the following non-factors to factors: condition, trial_within_condition, probabilitys191_p <- ggplot(s1_tr_summs_lc %>% filter(probability=="90-10"),
aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) +
geom_hline(yintercept = c(seq(.5, 1, .1))) +
geom_hline(yintercept = .5, size=3) +
geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") +
geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") +
geom_line(size=1, alpha=1) +
scale_color_manual(values=c("purple", "orange")) + ga + ap +
xlab("") + ylab("Proportion correct") + tol +
ggtitle("90-10") + tp + theme(axis.text.x = element_blank(), axis.ticks.x=element_blank()) +
scale_y_continuous(breaks=c(.6, .8, 1)) + ylim(.48, 1.02)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.s141_p <-ggplot(s1_tr_summs_lc %>% filter(probability=="40-10"),
aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) +
geom_hline(yintercept = c(seq(.5, 1, .1))) +
geom_hline(yintercept = .5, size=3) +
geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") +
geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") +
geom_line(size=1, alpha=1) +
scale_color_manual(values=c("purple", "orange")) + ga + ap +
xlab("Stim iteration") + ylab("") + tol +
ggtitle("40-10") + tp +
scale_y_continuous(breaks=c(.6, .8, 1)) + ylim(.48, 1.02)## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.s291_p <-ggplot(s2_tr_summs_lc %>% filter(probability=="90-10"),
aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) +
geom_hline(yintercept = c(seq(.5, 1, .1))) +
geom_hline(yintercept = .5, size=3) +
geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") +
geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") +
geom_line(size=1, alpha=1) +
scale_color_manual(values=c("purple", "orange")) + ga + ap +
xlab("") + ylab("") + tol +
ggtitle("90-10") + tp + theme(axis.text = element_blank(), axis.ticks=element_blank()) +
ylim(.48, 1.02)s241_p <-ggplot(s2_tr_summs_lc %>% filter(probability=="40-10"),
aes(x=as.numeric(trial_within_condition), y=m, group=condition, color=condition)) +
geom_hline(yintercept = c(seq(.5, 1, .1))) +
geom_hline(yintercept = .5, size=3) +
geom_vline(xintercept = 20, size=1, color="gray57", linetype="dotted") +
geom_ribbon(aes(x=as.numeric(trial_within_condition), ymin=m-se, ymax=m+se), alpha=.3, color="gray57") +
geom_line(size=1, alpha=1) +
scale_color_manual(values=c("purple", "orange")) + ga + ap +
xlab("Stim iteration") + ylab("") + tol +
ggtitle("40-10") + tp + theme(axis.text.y = element_blank(), axis.ticks.y=element_blank()) +
ylim(.48, 1.02)s1_lcs <- s191_p / s141_ps1_lcs## Warning: Removed 1 row containing missing values (`geom_line()`).
s2_lcs <- s291_p / s241_ps2_lcs
# ggsave("../paper/figs/fig_parts/s1_lcs.png", s1_lcs, height=6, width=5, dpi=700)
# ggsave("../paper/figs/fig_parts/s2_lcs.png", s2_lcs, height=6, width=5, dpi=700)Statistical analyses
Create sum contrast codes
# Specify overt as baseline factor so that negative effect corresponds to worst performance
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
#head(s1_train$cond_cc)
s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)
s1_sit$cond_cc <- factor(s1_sit$condition, levels=c("overt", "cognitive"))
contrasts(s1_sit$cond_cc) <- c(-.5, .5)
s2_sit$cond_cc <- factor(s2_sit$condition, levels=c("overt", "cognitive"))
contrasts(s2_sit$cond_cc) <- c(-.5, .5)Regressions of condition differences in performance
Experiment 1 — Learning and test models
summary(m1_train_full_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +
## scale(trial_within_condition) | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45054.2 45131.6 -22518.1 45036.2 39991
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.9157 -1.0686 0.4573 0.7047 1.2827
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4466 0.6683
## cond_cc1 0.4926 0.7018 -0.19
## scale(trial_within_condition) 0.1155 0.3399 0.82 0.01
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.04191 0.06123 17.017 < 2e-16 ***
## cond_cc1 -0.18054 0.06744 -2.677 0.00743 **
## scale(trial_within_condition) 0.36991 0.03300 11.209 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.174
## scl(trl_w_) 0.758 0.009summary(m1_test_full_regr <- glmer(correct ~ cond_cc + (cond_cc|ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7854.4 7889.3 -3922.2 7844.4 7995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1729 0.1319 0.2789 0.5791 1.2316
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.558 1.248
## cond_cc1 1.972 1.404 -0.17
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5690 0.1193 13.150 < 2e-16 ***
## cond_cc1 -0.4569 0.1510 -3.027 0.00247 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.173Experiment 2 — Learning and test models
summary(m2_train_full_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +
## scale(trial_within_condition) | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 48366.7 48445.0 -24174.4 48348.7 44151
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0444 -1.0336 0.4405 0.6872 1.4792
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6676 0.8171
## cond_cc1 0.5648 0.7515 -0.07
## scale(trial_within_condition) 0.1676 0.4094 0.87 -0.04
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.14912 0.07099 16.187 <2e-16 ***
## cond_cc1 -0.17465 0.06851 -2.549 0.0108 *
## scale(trial_within_condition) 0.39250 0.03732 10.517 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.071
## scl(trl_w_) 0.823 -0.034summary(m2_test_full_regr <- glmer(correct ~ cond_cc + (cond_cc|ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8426.3 8461.7 -4208.1 8416.3 8827
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2283 0.1250 0.2301 0.5751 1.5541
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.850 1.360
## cond_cc1 2.467 1.571 -0.15
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6720 0.1235 13.543 <2e-16 ***
## cond_cc1 -0.3543 0.1586 -2.234 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.156car::vif(m1_train_full_regr)## cond_cc scale(trial_within_condition)
## 1.000073 1.000073car::vif(m2_train_full_regr)## cond_cc scale(trial_within_condition)
## 1.001173 1.001173Robustness check 1: Same models with just random intercepts instead of both random intercepts and slopes
Experiment 1
summary(m1_train_ri_only_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(1 |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46071.1 46105.5 -23031.6 46063.1 39996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6755 -1.1011 0.5027 0.6818 1.2499
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.322 0.5675
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95856 0.05213 18.388 < 2e-16 ***
## cond_cc1 -0.14312 0.02271 -6.302 2.93e-10 ***
## scale(trial_within_condition) 0.28711 0.01147 25.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.006
## scl(trl_w_) 0.026 -0.008summary(m1_test_ri_only_regr <- glmer(correct ~ cond_cc + (1 |ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8080.6 8101.5 -4037.3 8074.6 7997
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4964 0.1516 0.3789 0.6174 1.2182
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.299 1.14
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.40802 0.10812 13.023 < 2e-16 ***
## cond_cc1 -0.31767 0.05567 -5.707 1.15e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.022Experiment 2
summary(m2_train_ri_only_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(1 | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 49635.6 49670.3 -24813.8 49627.6 44156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4739 -1.0679 0.4853 0.6668 1.2916
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4132 0.6428
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03285 0.05597 18.452 < 2e-16 ***
## cond_cc1 -0.15224 0.02195 -6.937 4.02e-12 ***
## scale(trial_within_condition) 0.27618 0.01108 24.921 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.007
## scl(trl_w_) 0.024 -0.008summary(m2_test_ri_only_regr <- glmer(correct ~ cond_cc + (1 |ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8769.8 8791.0 -4381.9 8763.8 8829
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6214 0.1441 0.3650 0.6075 1.2670
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.577 1.256
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.48452 0.11315 13.120 < 2e-16 ***
## cond_cc1 -0.20840 0.05335 -3.906 9.37e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.013Robustness check 2: Learning models with the time covariate removed
Experiments 1 and 2
summary(m1_train_no_cov_regr <- glmer(correct ~
cond_cc +
(cond_cc |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46182.9 46225.9 -23086.5 46172.9 39995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3753 -1.1134 0.5036 0.6746 1.0940
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.3342 0.5781
## cond_cc1 0.4616 0.6794 -0.20
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.97023 0.05308 18.279 < 2e-16 ***
## cond_cc1 -0.17612 0.06533 -2.696 0.00702 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.189summary(m2_train_no_cov_regr <- glmer(correct ~
cond_cc +
(cond_cc | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 49610.0 49653.5 -24800.0 49600.0 44155
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9719 -1.0789 0.4729 0.6802 1.0720
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4287 0.6548
## cond_cc1 0.5173 0.7192 -0.09
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0489 0.0570 18.401 < 2e-16 ***
## cond_cc1 -0.1697 0.0657 -2.583 0.00981 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.088sjPlot::tab_model(m1_train_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.83 | 2.51 – 3.20 | <0.001 |
| cond cc1 | 0.83 | 0.73 – 0.95 | 0.007 |
| trial within condition | 1.45 | 1.36 – 1.54 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.45 | ||
| τ11 ID.cond_cc1 | 0.49 | ||
| τ11 ID.scale(trial_within_condition) | 0.12 | ||
| ρ01 | -0.19 | ||
| 0.82 | |||
| ICC | 0.17 | ||
| N ID | 125 | ||
| Observations | 40000 | ||
| Marginal R2 / Conditional R2 | 0.035 / 0.202 | ||
sjPlot::tab_model(m2_train_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 3.16 | 2.75 – 3.63 | <0.001 |
| cond cc1 | 0.84 | 0.73 – 0.96 | 0.011 |
| trial within condition | 1.48 | 1.38 – 1.59 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.67 | ||
| τ11 ID.cond_cc1 | 0.56 | ||
| τ11 ID.scale(trial_within_condition) | 0.17 | ||
| ρ01 | -0.07 | ||
| 0.87 | |||
| ICC | 0.23 | ||
| N ID | 138 | ||
| Observations | 44160 | ||
| Marginal R2 / Conditional R2 | 0.037 / 0.257 | ||
sjPlot::tab_model(m1_test_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 4.80 | 3.80 – 6.07 | <0.001 |
| cond cc1 | 0.63 | 0.47 – 0.85 | 0.002 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.56 | ||
| τ11 ID.cond_cc1 | 1.97 | ||
| ρ01 ID | -0.17 | ||
| ICC | 0.38 | ||
| N ID | 125 | ||
| Observations | 8000 | ||
| Marginal R2 / Conditional R2 | 0.010 / 0.390 | ||
sjPlot::tab_model(m2_test_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 4.18 – 6.78 | <0.001 |
| cond cc1 | 0.70 | 0.51 – 0.96 | 0.025 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.85 | ||
| τ11 ID.cond_cc1 | 2.47 | ||
| ρ01 ID | -0.15 | ||
| ICC | 0.43 | ||
| N ID | 138 | ||
| Observations | 8832 | ||
| Marginal R2 / Conditional R2 | 0.005 / 0.432 | ||
Robustness check 3: Covariate for answer type
s1_train[s1_train$response=="sum", "resp_type_cc_cog"] <- .5
s1_train[s1_train$response=="difference", "resp_type_cc_cog"] <- -.5
s1_train[s1_train$response=="top", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="bottom", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="top", "resp_type_cc_overt"] <- .5
s1_train[s1_train$response=="bottom", "resp_type_cc_overt"] <- -.5
s1_train[s1_train$response=="sum", "resp_type_cc_overt"] <- 0
s1_train[s1_train$response=="difference", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_cog"] <- .5
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_cog"] <- -.5
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_overt"] <- .5
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_overt"] <- -.5
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="alphabetize", "resp_type_cc_cog"] <- .5
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_train[s2_train$response=="slash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="backslash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="slash", "resp_type_cc_overt"] <- .5
s2_train[s2_train$response=="backslash", "resp_type_cc_overt"] <- -.5
s2_train[s2_train$response=="alphabetize", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_cog"] <- .5
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_overt"] <- .5
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_overt"] <- -.5
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_overt"] <- 0
table(s1_train$resp_type_cc_cog)##
## -0.5 0 0.5
## 9020 20000 10980table(s1_train$resp_type_cc_overt)##
## -0.5 0 0.5
## 9451 20000 10549table(s1_sit$resp_type_cc_cog)##
## -0.5 0 0.5
## 1851 4000 2149table(s1_sit$resp_type_cc_overt)##
## -0.5 0 0.5
## 1975 4000 2025table(s2_train$resp_type_cc_cog)##
## -0.5 0 0.5
## 10182 22080 11898table(s2_train$resp_type_cc_overt)##
## -0.5 0 0.5
## 11186 22080 10894table(s2_sit$resp_type_cc_cog)##
## -0.5 0 0.5
## 2024 4416 2392table(s2_sit$resp_type_cc_overt)##
## -0.5 0 0.5
## 2239 4416 2177Separately coded contrast variables for cog and overt. Random slopes of these removed for test bc singular (which makes sense bc of few observations/subj)
summary(m1_train_full_regr_at <- glmer(correct ~
cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +
## resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +
## resp_type_cc_cog + resp_type_cc_overt | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45040 45212 -22500 45000 39980
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2192 -1.0643 0.4553 0.7031 1.4380
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.44641 0.6681
## cond_cc1 0.49626 0.7045 -0.19
## scale(trial_within_condition) 0.11598 0.3406 0.81 0.01
## resp_type_cc_cog 0.05453 0.2335 0.19 -0.31 0.33
## resp_type_cc_overt 0.07932 0.2816 -0.04 0.23 0.16 0.70
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.06182 0.06142 17.287 < 2e-16 ***
## cond_cc1 -0.18292 0.06840 -2.674 0.00749 **
## scale(trial_within_condition) 0.37205 0.03308 11.247 < 2e-16 ***
## resp_type_cc_cog -0.07227 0.04151 -1.741 0.08165 .
## resp_type_cc_overt -0.08410 0.04532 -1.855 0.06353 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1 -0.172
## scl(trl_w_) 0.753 0.004
## rsp_typ_cc_c 0.077 -0.175 0.164
## rsp_typ_cc_v -0.034 0.137 0.088 0.197# Singular
# summary(m1_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
# (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID),
# data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m1_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc |ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7857.3 7906.2 -3921.6 7843.3 7993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0858 0.1296 0.2776 0.5823 1.2399
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.555 1.247
## cond_cc1 1.979 1.407 -0.17
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.57061 0.11923 13.173 < 2e-16 ***
## cond_cc1 -0.45590 0.15119 -3.015 0.00257 **
## resp_type_cc_cog -0.05969 0.08348 -0.715 0.47460
## resp_type_cc_overt -0.07056 0.08935 -0.790 0.42972
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1 -0.175
## rsp_typ_cc_c -0.015 -0.021
## rsp_typ_cc_v -0.007 0.011 0.000summary(m2_train_full_regr_at <- glmer(correct ~
cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +
## resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +
## resp_type_cc_cog + resp_type_cc_overt | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 48297.6 48471.5 -24128.8 48257.6 44140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.4618 -1.0313 0.4323 0.6859 1.4685
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6756 0.8220
## cond_cc1 0.5985 0.7736 -0.07
## scale(trial_within_condition) 0.1666 0.4081 0.87 -0.04
## resp_type_cc_cog 0.2001 0.4474 0.01 -0.26 0.01
## resp_type_cc_overt 0.1730 0.4160 -0.27 0.00 -0.27 0.18
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.194005 0.071628 16.670 <2e-16 ***
## cond_cc1 -0.163191 0.071098 -2.295 0.0217 *
## scale(trial_within_condition) 0.398418 0.037258 10.693 <2e-16 ***
## resp_type_cc_cog -0.122259 0.052879 -2.312 0.0208 *
## resp_type_cc_overt -0.001977 0.051645 -0.038 0.9695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1 -0.064
## scl(trl_w_) 0.821 -0.037
## rsp_typ_cc_c -0.001 -0.199 0.016
## rsp_typ_cc_v -0.183 0.007 -0.170 0.087# Singular
# summary(m2_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
# (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID),
# data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m2_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc |ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8427.9 8477.5 -4207.0 8413.9 8825
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3845 0.1211 0.2303 0.5752 1.5798
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.848 1.359
## cond_cc1 2.469 1.571 -0.16
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67413 0.12343 13.564 <2e-16 ***
## cond_cc1 -0.35472 0.15863 -2.236 0.0253 *
## resp_type_cc_cog -0.05734 0.08376 -0.685 0.4936
## resp_type_cc_overt 0.11865 0.08580 1.383 0.1667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1 -0.159
## rsp_typ_cc_c -0.014 -0.022
## rsp_typ_cc_v 0.010 -0.015 0.000Robustness check 4: paired t-test on arcsine-transformed proportion data instead of mixed-effects model
s1_summs_id <- s1_train %>% group_by(condition, ID) %>%
summarize(m=mean(correct)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_summs_id <- s2_train %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_summs_test_id <- s1_sit %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_summs_test_id <- s2_sit %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.t.test(asin(sqrt(m)) ~ condition, data=s1_summs_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.5152, df = 124, p-value = 0.01318
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.061942537 -0.007385729
## sample estimates:
## mean difference
## -0.03466413t.test(asin(sqrt(m)) ~ condition, data=s2_summs_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.5139, df = 137, p-value = 0.0131
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.061281292 -0.007319414
## sample estimates:
## mean difference
## -0.03430035t.test(asin(sqrt(m)) ~ condition, data=s1_summs_test_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.6957, df = 124, p-value = 0.008001
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.12657517 -0.01939599
## sample estimates:
## mean difference
## -0.07298558t.test(asin(sqrt(m)) ~ condition, data=s2_summs_test_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -1.8641, df = 137, p-value = 0.06445
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.105733994 0.003119297
## sample estimates:
## mean difference
## -0.05130735Additional check for resub: Is there moderation by difference in reaction times?
mean(s1_train[s1_train$condition=="cognitive", "rt_key_1"])-mean(s1_train[s1_train$condition=="overt", "rt_key_1"])## [1] 246.5527mean(s1_sit[s1_sit$condition=="cognitive", "rt_key_1"])-mean(s1_sit[s1_sit$condition=="overt", "rt_key_1"])## [1] 188.8504mean(s2_train[s2_train$condition=="cognitive", "rt_key_1"])-mean(s2_train[s2_train$condition=="overt", "rt_key_1"])## [1] 680.3651mean(s2_sit[s2_sit$condition=="cognitive", "rt_key_1"])-mean(s2_sit[s2_sit$condition=="overt", "rt_key_1"])## [1] 643.3419s1_tr_rt_inds <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_tr_diffs_inds <-
data.frame("diff"=s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "m"]-s1_tr_rt_inds[s1_tr_rt_inds$condition=="overt", "m"],
"ID"=s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "ID"])
assert(all(s1_tr_rt_inds[s1_tr_rt_inds$condition=="cognitive", "ID"]==s1_tr_rt_inds[s1_tr_rt_inds$condition=="overt", "ID"]))
s1_te_rt_inds <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_te_diffs_inds <-
data.frame("diff"=s1_te_rt_inds[s1_te_rt_inds$condition=="cognitive", "m"]-s1_te_rt_inds[s1_te_rt_inds$condition=="overt", "m"],
"ID"=s1_te_rt_inds[s1_te_rt_inds$condition=="cognitive", "ID"])
s2_tr_rt_inds <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_tr_diffs_inds <-
data.frame("diff"=s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "m"]-s2_tr_rt_inds[s2_tr_rt_inds$condition=="overt", "m"],
"ID"=s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "ID"])
assert(all(s2_tr_rt_inds[s2_tr_rt_inds$condition=="cognitive", "ID"]==s2_tr_rt_inds[s2_tr_rt_inds$condition=="overt", "ID"]))
s2_te_rt_inds <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(rt_key_1)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_te_diffs_inds <-
data.frame("diff"=s2_te_rt_inds[s2_te_rt_inds$condition=="cognitive", "m"]-s2_te_rt_inds[s2_te_rt_inds$condition=="overt", "m"],
"ID"=s2_te_rt_inds[s2_te_rt_inds$condition=="cognitive", "ID"])Add the RT diffs to dataframes
s1_IDs <- unique(s1_train$ID)
for (i in 1:length(s1_IDs)) {
s1_train[s1_train$ID == s1_IDs[i], "rt_diff"] <- s1_tr_diffs_inds[s1_tr_diffs_inds$ID == s1_IDs[i], "m"]
s1_sit[s1_sit$ID == s1_IDs[i], "rt_diff"] <- s1_te_diffs_inds[s1_te_diffs_inds$ID == s1_IDs[i], "m"]
}
# Spot checks
# s1_train %>% filter(ID==25) %>% select(rt_diff)
# s1_tr_diffs_inds %>% filter(ID==25) %>% select(m)
#
# s1_sit %>% filter(ID==25) %>% select(rt_diff)
# s1_te_diffs_inds %>% filter(ID==25) %>% select(m)
#
# s1_train %>% filter(ID==102) %>% select(rt_diff)
# s1_tr_diffs_inds %>% filter(ID==102) %>% select(m)
#
# s1_sit %>% filter(ID==102) %>% select(rt_diff)
# s1_te_diffs_inds %>% filter(ID==102) %>% select(m)s2_IDs <- unique(s2_train$ID)
for (i in 1:length(s1_IDs)) {
s2_train[s2_train$ID == s2_IDs[i], "rt_diff"] <- s2_tr_diffs_inds[s2_tr_diffs_inds$ID == s2_IDs[i], "m"]
s2_sit[s2_sit$ID == s2_IDs[i], "rt_diff"] <- s2_te_diffs_inds[s2_te_diffs_inds$ID == s2_IDs[i], "m"]
}
# Spot checks
# s2_train %>% filter(ID==25) %>% select(rt_diff)
# s2_tr_diffs_inds %>% filter(ID==25) %>% select(m)
#
# s2_sit %>% filter(ID==25) %>% select(rt_diff)
# s2_te_diffs_inds %>% filter(ID==25) %>% select(m)
#
# s2_train %>% filter(ID==102) %>% select(rt_diff)
# s2_tr_diffs_inds %>% filter(ID==102) %>% select(m)
# #
# s2_sit %>% filter(ID==102) %>% select(rt_diff)
# s2_te_diffs_inds %>% filter(ID==102) %>% select(m)Examine as potential moderator
Reduced complexity where needed because more complex models were singular
summary(m1_train_full_regr <- glmer(correct ~
cond_cc*scale(rt_diff) + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + scale(trial_within_condition) +
## (cond_cc + scale(trial_within_condition) | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45057.4 45152.0 -22517.7 45035.4 39989
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.8320 -1.0683 0.4574 0.7050 1.2730
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4451 0.6672
## cond_cc1 0.4900 0.7000 -0.18
## scale(trial_within_condition) 0.1156 0.3400 0.82 0.01
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.04184 0.06112 17.046 < 2e-16 ***
## cond_cc1 -0.18034 0.06728 -2.681 0.00735 **
## scale(rt_diff) 0.02641 0.03960 0.667 0.50490
## scale(trial_within_condition) 0.36992 0.03300 11.209 < 2e-16 ***
## cond_cc1:scale(rt_diff) -0.04805 0.06688 -0.718 0.47249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 scl(_) sc(__)
## cond_cc1 -0.170
## scl(rt_dff) 0.000 0.000
## scl(trl_w_) 0.759 0.010 0.000
## cnd_cc1:(_) 0.000 -0.001 -0.274 0.000summary(m1_test_full_regr <- glmer(correct ~ cond_cc*scale(rt_diff) + (cond_cc |ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + (cond_cc | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7856.2 7905.1 -3921.1 7842.2 7993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0633 0.1280 0.2738 0.5815 1.2273
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.537 1.240
## cond_cc1 1.945 1.395 -0.15
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5680 0.1186 13.221 < 2e-16 ***
## cond_cc1 -0.4534 0.1502 -3.018 0.00254 **
## scale(rt_diff) 0.1212 0.1159 1.046 0.29555
## cond_cc1:scale(rt_diff) -0.1729 0.1412 -1.225 0.22051
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 scl(_)
## cond_cc1 -0.160
## scl(rt_dff) 0.006 -0.015
## cnd_cc1:(_) -0.016 0.015 -0.145Reduced complexity for these two because more complex models singular
summary(m2_train_full_regr <- glmer(correct ~
cond_cc*scale(rt_diff) + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) |ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + scale(trial_within_condition) +
## (cond_cc + scale(trial_within_condition) | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 43489.9 43584.5 -21734.0 43467.9 39989
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.2304 -1.0286 0.4327 0.6848 1.4586
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.7250 0.8515
## cond_cc1 0.5118 0.7154 -0.08
## scale(trial_within_condition) 0.1803 0.4247 0.88 -0.07
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.17379 0.07769 15.109 < 2e-16 ***
## cond_cc1 -0.15927 0.06908 -2.306 0.02113 *
## scale(rt_diff) 0.05355 0.04297 1.246 0.21271
## scale(trial_within_condition) 0.40406 0.04059 9.954 < 2e-16 ***
## cond_cc1:scale(rt_diff) -0.22740 0.06887 -3.302 0.00096 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 scl(_) sc(__)
## cond_cc1 -0.074
## scl(rt_dff) 0.001 -0.007
## scl(trl_w_) 0.832 -0.064 -0.001
## cnd_cc1:(_) -0.006 0.006 -0.032 -0.002summary(m2_test_full_regr <- glmer(correct ~ cond_cc*scale(rt_diff) + (cond_cc |ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * scale(rt_diff) + (cond_cc | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7490.3 7539.2 -3738.2 7476.3 7993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6644 0.1168 0.2216 0.5651 1.5704
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.904 1.380
## cond_cc1 2.307 1.519 -0.09
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7299 0.1318 13.124 <2e-16 ***
## cond_cc1 -0.2647 0.1640 -1.614 0.1066
## scale(rt_diff) 0.2416 0.1304 1.852 0.0640 .
## cond_cc1:scale(rt_diff) -0.3419 0.1580 -2.164 0.0305 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 scl(_)
## cond_cc1 -0.093
## scl(rt_dff) 0.020 -0.018
## cnd_cc1:(_) -0.017 0.052 -0.100Additional check for resub 2: Is performance better when the first learning block of the condition happens second, reflective of meta-learning
s1_with_cond_order <- read.csv("../data/cleaned_files/s1_train_deident_with_first_cond.csv")
s2_with_cond_order <- read.csv("../data/cleaned_files/s2_train_with_exclusions_IDed_with_first_cond_wdeID.csv")Put the first assignments into regular dfs
unique_IDs <- unique(s1_train$ID)
for (s in 1:length(unique_IDs)) {
s1_train[s1_train$ID == unique_IDs[s], "first_cond"] <-
s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"]
s1_sit[s1_sit$ID == unique_IDs[s], "first_cond"] <-
unique(s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"])
}s1_train %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: condition [2]
## condition first_cond m
## <chr> <chr> <dbl>
## 1 cognitive cognTraining 0.670
## 2 cognitive overtTraining 0.715
## 3 overt cognTraining 0.723
## 4 overt overtTraining 0.716s1_sit %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: condition [2]
## condition first_cond m
## <chr> <chr> <dbl>
## 1 cognitive cognTraining 0.703
## 2 cognitive overtTraining 0.752
## 3 overt cognTraining 0.790
## 4 overt overtTraining 0.765unique_IDs <- unique(s2_train$ID)
for (s in 1:length(unique_IDs)) {
s2_train[s2_train$ID == unique_IDs[s], "first_cond"] <-
s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"]
s2_sit[s2_sit$ID == unique_IDs[s], "first_cond"] <-
unique(s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"])
}s2_train %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: condition [2]
## condition first_cond m
## <chr> <chr> <dbl>
## 1 cognitive cognitive 0.667
## 2 cognitive overt 0.733
## 3 overt cognitive 0.734
## 4 overt overt 0.728s2_sit %>% group_by(condition, first_cond) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.## # A tibble: 4 × 3
## # Groups: condition [2]
## condition first_cond m
## <chr> <chr> <dbl>
## 1 cognitive cognitive 0.687
## 2 cognitive overt 0.784
## 3 overt cognitive 0.775
## 4 overt overt 0.769Random assignment of first condition was largely even, with more assignments to cognitive first in study 1 and overt first in condition 2
table(s1_train$first_cond)##
## cognTraining overtTraining
## 20480 19520table(s2_train$first_cond)##
## cognitive overt
## 20480 23680# Specify overt as baseline factor so that negative effect corresponds to worst performance
s1_train$fc_cond <- factor(s1_train$first_cond, levels=c("overtTraining", "cognTraining"))
contrasts(s1_train$fc_cond) <- c(-.5, .5)
#head(s1_train$fc_cond)
s2_train$fc_cond <- factor(s2_train$first_cond, levels=c("overt", "cognitive"))
contrasts(s2_train$fc_cond) <- c(-.5, .5)
s1_sit$fc_cond <- factor(s1_sit$first_cond, levels=c("overtTraining", "cognTraining"))
contrasts(s1_sit$fc_cond) <- c(-.5, .5)
s2_sit$fc_cond <- factor(s2_sit$first_cond, levels=c("overt", "cognitive"))
contrasts(s2_sit$fc_cond) <- c(-.5, .5)Consistent pattern whereby the condition difference is most pronounced when cognitive is first; main effect remains significant
# Reduced RE from interaction bc of non-convergence with full
summary(m1_train_full_regr <- glmer(correct ~
cond_cc*fc_cond + scale(trial_within_condition) +
(cond_cc + fc_cond + scale(trial_within_condition) |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * fc_cond + scale(trial_within_condition) +
## (cond_cc + fc_cond + scale(trial_within_condition) | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45053.5 45182.5 -22511.8 45023.5 39985
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.9434 -1.0673 0.4574 0.7052 1.2873
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.3225 0.5679
## cond_cc1 0.4698 0.6854 -0.27
## fc_cond1 0.5118 0.7154 -0.30 0.00
## scale(trial_within_condition) 0.1154 0.3398 0.96 -0.02 -0.22
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03334 0.06137 16.839 < 2e-16 ***
## cond_cc1 -0.17725 0.06609 -2.682 0.00732 **
## fc_cond1 -0.02511 0.08103 -0.310 0.75662
## scale(trial_within_condition) 0.36981 0.03298 11.212 < 2e-16 ***
## cond_cc1:fc_cond1 -0.29105 0.13287 -2.191 0.02849 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 fc_cn1 sc(__)
## cond_cc1 -0.210
## fc_cond1 -0.285 0.006
## scl(trl_w_) 0.758 -0.015 -0.152
## cnd_cc1:f_1 0.001 -0.025 -0.301 -0.003# Reduced for same reason
summary(m1_test_full_regr <- glmer(correct ~ cond_cc*fc_cond + (cond_cc + fc_cond|ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * fc_cond + (cond_cc + fc_cond | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7847.1 7917.0 -3913.6 7827.1 7990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2749 0.1146 0.2742 0.5759 1.2486
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.197 1.094
## cond_cc1 2.004 1.416 -0.14
## fc_cond1 1.889 1.374 -0.23 -0.64
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5963 0.1240 12.874 < 2e-16 ***
## cond_cc1 -0.4114 0.1520 -2.707 0.00679 **
## fc_cond1 -0.1968 0.2456 -0.801 0.42295
## cond_cc1:fc_cond1 -0.9458 0.3070 -3.081 0.00207 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 fc_cn1
## cond_cc1 -0.113
## fc_cond1 -0.224 -0.299
## cnd_cc1:f_1 -0.314 -0.063 -0.106# Run with full interaction because
summary(m2_train_full_regr <- glmer(correct ~
cond_cc*fc_cond + scale(trial_within_condition) +
(cond_cc*fc_cond + scale(trial_within_condition) | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * fc_cond + scale(trial_within_condition) +
## (cond_cc * fc_cond + scale(trial_within_condition) | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 48362.6 48536.6 -24161.3 48322.6 44140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.2482 -1.0364 0.4408 0.6878 1.4992
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.5537 0.7441
## cond_cc1 0.3106 0.5573 -0.28
## fc_cond1 0.4216 0.6493 -0.11 -0.30
## scale(trial_within_condition) 0.1689 0.4110 0.94 -0.19 -0.16
## cond_cc1:fc_cond1 0.8624 0.9287 -0.29 -0.10 0.14 -0.23
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.14228 0.07120 16.042 < 2e-16 ***
## cond_cc1 -0.21284 0.06739 -3.158 0.00159 **
## fc_cond1 0.02138 0.08324 0.257 0.79734
## scale(trial_within_condition) 0.39290 0.03745 10.491 < 2e-16 ***
## cond_cc1:fc_cond1 -0.52759 0.13473 -3.916 9.01e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 fc_cn1 sc(__)
## cond_cc1 -0.179
## fc_cond1 0.007 -0.254
## scl(trl_w_) 0.813 -0.133 -0.088
## cnd_cc1:f_1 -0.262 -0.001 -0.109 -0.132summary(m2_test_full_regr <- glmer(correct ~ cond_cc*fc_cond + (cond_cc*fc_cond|ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc * fc_cond + (cond_cc * fc_cond | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8430.5 8529.7 -4201.3 8402.5 8818
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2679 0.1224 0.2404 0.5711 1.5615
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.411 1.188
## cond_cc1 1.957 1.399 -0.21
## fc_cond1 1.651 1.285 -0.28 -0.07
## cond_cc1:fc_cond1 1.458 1.208 -0.13 -0.05 -0.10
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6640 0.1221 13.625 < 2e-16 ***
## cond_cc1 -0.3674 0.1554 -2.364 0.01808 *
## fc_cond1 -0.4580 0.2441 -1.876 0.06068 .
## cond_cc1:fc_cond1 -0.8071 0.3107 -2.598 0.00939 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 fc_cn1
## cond_cc1 -0.183
## fc_cond1 -0.158 -0.098
## cnd_cc1:f_1 -0.098 -0.007 -0.183For resub, correlation between training and test performance
cor.test(
as.numeric(unlist(s1_train_summs[s1_train_summs$condition=="overt", "m"])),
as.numeric(unlist(s1_test_summs[s1_test_summs$condition=="overt", "m"]))
)##
## Pearson's product-moment correlation
##
## data: as.numeric(unlist(s1_train_summs[s1_train_summs$condition == "overt", "m"])) and as.numeric(unlist(s1_test_summs[s1_test_summs$condition == "overt", "m"]))
## t = 12.528, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6599115 0.8169410
## sample estimates:
## cor
## 0.7487498cor.test(
as.numeric(unlist(s1_train_summs[s1_train_summs$condition=="cognitive", "m"])),
as.numeric(unlist(s1_test_summs[s1_test_summs$condition=="cognitive", "m"]))
)##
## Pearson's product-moment correlation
##
## data: as.numeric(unlist(s1_train_summs[s1_train_summs$condition == "cognitive", "m"])) and as.numeric(unlist(s1_test_summs[s1_test_summs$condition == "cognitive", "m"]))
## t = 14.454, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7177575 0.8504835
## sample estimates:
## cor
## 0.7933662cor.test(
as.numeric(unlist(s2_train_summs[s2_train_summs$condition=="overt", "m"])),
as.numeric(unlist(s2_test_summs[s2_test_summs$condition=="overt", "m"]))
)##
## Pearson's product-moment correlation
##
## data: as.numeric(unlist(s2_train_summs[s2_train_summs$condition == "overt", "m"])) and as.numeric(unlist(s2_test_summs[s2_test_summs$condition == "overt", "m"]))
## t = 11.385, df = 136, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6016955 0.7751783
## sample estimates:
## cor
## 0.6985625cor.test(
as.numeric(unlist(s2_train_summs[s2_train_summs$condition=="cognitive", "m"])),
as.numeric(unlist(s2_test_summs[s2_test_summs$condition=="cognitive", "m"]))
)##
## Pearson's product-moment correlation
##
## data: as.numeric(unlist(s2_train_summs[s2_train_summs$condition == "cognitive", "m"])) and as.numeric(unlist(s2_test_summs[s2_test_summs$condition == "cognitive", "m"]))
## t = 13.631, df = 136, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6789692 0.8225245
## sample estimates:
## cor
## 0.7598609Delay effects
Create sum contrast codes for delay
# Specify overt as baseline factor so that negative effect corresponds to worst performance
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)#head(factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay")))
s1_train$delay_cc <- factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s1_train$delay_cc) <- c(-.5, .5)
s2_train$delay_cc <- factor(if_else(s2_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s2_train$delay_cc) <- c(-.5, .5)Building up in complexity…
Effect of delay with no moderation
summary(m1_delay <- glmer(correct ~ delay_cc + (delay_cc|ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 43647.9 43690.6 -21818.9 43637.9 37995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2305 -1.1242 0.5124 0.6818 1.0042
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.32784 0.5726
## delay_cc1 0.02447 0.1564 0.43
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.07222 0.05324 20.14 <2e-16 ***
## delay_cc1 -0.34260 0.03240 -10.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## delay_cc1 0.043summary(m2_delay <- glmer(correct ~ delay_cc + (delay_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 47084.2 47127.5 -23537.1 47074.2 41947
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2244 -1.1013 0.4908 0.6573 1.0211
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.42730 0.6537
## delay_cc1 0.01198 0.1095 0.44
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.13294 0.05749 19.71 <2e-16 ***
## delay_cc1 -0.30347 0.03024 -10.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## delay_cc1 -0.001Now moderated by condition
#Singular for both
# summary(m1_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc*cond_cc|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# summary(m2_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc*cond_cc|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Simplified
summary(m1_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc + cond_cc|ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 43090.3 43175.7 -21535.1 43070.3 37990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0844 -1.0664 0.4738 0.6770 1.1528
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.35727 0.5977
## delay_cc1 0.02429 0.1559 0.37
## cond_cc1 0.50501 0.7106 -0.22 -0.17
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.10896 0.05554 19.967 <2e-16 ***
## delay_cc1 -0.35536 0.03271 -10.864 <2e-16 ***
## cond_cc1 -0.17949 0.07027 -2.554 0.0106 *
## delay_cc1:cond_cc1 -0.01904 0.05747 -0.331 0.7404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1
## delay_cc1 0.019
## cond_cc1 -0.196 -0.055
## dly_cc1:c_1 0.006 -0.041 -0.226summary(m2_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc + cond_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46388.5 46474.9 -23184.2 46368.5 41942
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.9099 -1.0345 0.4587 0.6672 1.1512
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.45639 0.6756
## delay_cc1 0.01049 0.1024 0.55
## cond_cc1 0.56015 0.7484 -0.09 -0.21
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.17084 0.05936 19.725 <2e-16 ***
## delay_cc1 -0.30464 0.03032 -10.046 <2e-16 ***
## cond_cc1 -0.13973 0.07005 -1.995 0.0461 *
## delay_cc1:cond_cc1 -0.13157 0.05525 -2.381 0.0172 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1
## delay_cc1 0.023
## cond_cc1 -0.085 -0.054
## dly_cc1:c_1 0.001 -0.025 -0.214Model including appropriate covariates at the highest complexity that would fit (removed trial-within-cond random slope and delay*cond slope due to singular errors in studies 2 and 1 respectively)
# summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc + cond_cc + scale(trial_within_condition)|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# # Singular
# summary(m2_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc + cond_cc + scale(trial_within_condition)|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Also singular
# summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc*cond_cc + 1|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
(delay_cc + cond_cc |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +
## (delay_cc + cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 42576.5 42670.5 -21277.2 42554.5 37989
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2409 -1.0187 0.4627 0.6726 1.4465
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.36544 0.6045
## delay_cc1 0.02941 0.1715 0.36
## cond_cc1 0.51948 0.7208 -0.21 -0.18
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.11508 0.05619 19.845 < 2e-16 ***
## delay_cc1 -0.36388 0.03352 -10.857 < 2e-16 ***
## cond_cc1 -0.18364 0.07121 -2.579 0.00991 **
## scale(trial_within_condition) 0.27402 0.01215 22.553 < 2e-16 ***
## delay_cc1:cond_cc1 -0.01022 0.05790 -0.177 0.85983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1 0.029
## cond_cc1 -0.194 -0.064
## scl(trl_w_) 0.018 -0.019 -0.004
## dly_cc1:c_1 0.007 -0.041 -0.225 0.006summary(m2_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
(delay_cc + cond_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +
## (delay_cc + cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45878.5 45973.6 -22928.2 45856.5 41941
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8197 -0.9779 0.4437 0.6526 1.4299
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.46804 0.6841
## delay_cc1 0.01325 0.1151 0.47
## cond_cc1 0.57474 0.7581 -0.09 -0.18
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.17603 0.06016 19.549 <2e-16 ***
## delay_cc1 -0.30940 0.03083 -10.036 <2e-16 ***
## cond_cc1 -0.14542 0.07100 -2.048 0.0405 *
## scale(trial_within_condition) 0.26368 0.01172 22.489 <2e-16 ***
## delay_cc1:cond_cc1 -0.11654 0.05571 -2.092 0.0364 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1 0.016
## cond_cc1 -0.084 -0.052
## scl(trl_w_) 0.016 -0.013 -0.004
## dly_cc1:c_1 0.001 -0.025 -0.214 0.010car::vif(m1_full_delay_cond)## delay_cc cond_cc
## 1.007763 1.059556
## scale(trial_within_condition) delay_cc:cond_cc
## 1.000389 1.056908car::vif(m2_full_delay_cond)## delay_cc cond_cc
## 1.004255 1.051818
## scale(trial_within_condition) delay_cc:cond_cc
## 1.000268 1.049647For resub, rerun interaction with continuous delay
summary(m1_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
(cond_cc|ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(delay) * cond_cc + scale(trial_within_condition) +
## (cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 42668.7 42737.1 -21326.4 42652.7 37992
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4888 -1.0349 0.4694 0.6753 1.5549
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.3873 0.6224
## cond_cc1 0.5156 0.7181 -0.22
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.02230 0.05713 17.894 < 2e-16 ***
## scale(delay) -0.11266 0.01154 -9.766 < 2e-16 ***
## cond_cc1 -0.18989 0.06914 -2.746 0.00603 **
## scale(trial_within_condition) 0.27169 0.01212 22.421 < 2e-16 ***
## scale(delay):cond_cc1 0.02219 0.02307 0.962 0.33610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(d) cnd_c1 sc(__)
## scale(dely) -0.009
## cond_cc1 -0.208 0.003
## scl(trl_w_) 0.015 -0.015 -0.003
## scl(dly):_1 0.002 -0.018 -0.015 -0.001# Singular
# summary(m1_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
# (scale(delay) + cond_cc|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))# Singular
# summary(m1_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
# (scale(delay)*cond_cc |ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))summary(m2_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
( cond_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ scale(delay) * cond_cc + scale(trial_within_condition) +
## (cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45932.2 46001.3 -22958.1 45916.2 41944
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8928 -0.9866 0.4450 0.6548 1.5689
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4873 0.6981
## cond_cc1 0.5749 0.7582 -0.10
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.097613 0.060775 18.060 <2e-16 ***
## scale(delay) -0.108518 0.011140 -9.741 <2e-16 ***
## cond_cc1 -0.177724 0.069359 -2.562 0.0104 *
## scale(trial_within_condition) 0.263420 0.011711 22.494 <2e-16 ***
## scale(delay):cond_cc1 -0.004329 0.022280 -0.194 0.8459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) scl(d) cnd_c1 sc(__)
## scale(dely) -0.009
## cond_cc1 -0.092 0.001
## scl(trl_w_) 0.014 -0.010 -0.003
## scl(dly):_1 0.001 -0.030 -0.015 -0.003# Also singular
# summary(m2_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
# (scale(delay) + cond_cc|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# # singular
# summary(m2_full_delay_cond_cont <- glmer(correct ~ scale(delay)*cond_cc + scale(trial_within_condition) +
# (scale(delay)*cond_cc|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))Heterogeneity in condition effect
Calculate some bxal differences
Study 1
s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
s1_perf[s1_perf$condition=="cognitive", "m"]))
# .. and test subject-level differences in performance
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))
s1_perf_overall <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf_overall <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))
# Modest correlations across conditions
cor.test(unlist(s1_perf[s1_perf$condition=="overt","m"]), unlist(s1_perf[s1_perf$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s1_perf[s1_perf$condition == "overt", "m"]) and unlist(s1_perf[s1_perf$condition == "cognitive", "m"])
## t = 5.461, df = 123, p-value = 2.507e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2885290 0.5729171
## sample estimates:
## cor
## 0.4417538cor.test(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]), unlist(s1_perf_test[s1_perf_test$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s1_perf_test[s1_perf_test$condition == "overt", "m"]) and unlist(s1_perf_test[s1_perf_test$condition == "cognitive", "m"])
## t = 5.6112, df = 123, p-value = 1.263e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2995928 0.5809966
## sample estimates:
## cor
## 0.4514492Examining if between-condition diffs are more common at the highest level of perf, because at lower ends more noise in bx. But not much relationship
cor.test(s1_perf_overall$m, o_min_c_perf_s1)##
## Pearson's product-moment correlation
##
## data: s1_perf_overall$m and o_min_c_perf_s1
## t = 1.0074, df = 123, p-value = 0.3157
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.08652362 0.26190545
## sample estimates:
## cor
## 0.09045834ComparePars(s1_perf_overall$m, o_min_c_perf_s1, use_identity_line = 0)
Study 2
s2_perf <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s2 <- as.numeric(unlist(s2_perf[s2_perf$condition=="overt","m"]-
s2_perf[s2_perf$condition=="cognitive", "m"]))
s2_perf_test <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]-
s2_perf_test[s2_perf_test$condition=="cognitive", "m"]))
s2_perf_overall <- s2_train %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_overall <- s2_sit %>% group_by(ID) %>% summarize(m=mean(correct))
cor.test(unlist(s2_perf[s2_perf$condition=="overt","m"]), unlist(s2_perf[s2_perf$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s2_perf[s2_perf$condition == "overt", "m"]) and unlist(s2_perf[s2_perf$condition == "cognitive", "m"])
## t = 6.0926, df = 136, p-value = 1.074e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3207662 0.5848974
## sample estimates:
## cor
## 0.4630508cor.test(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]), unlist(s2_perf_test[s2_perf_test$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s2_perf_test[s2_perf_test$condition == "overt", "m"]) and unlist(s2_perf_test[s2_perf_test$condition == "cognitive", "m"])
## t = 4.8527, df = 136, p-value = 3.29e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2319669 0.5180282
## sample estimates:
## cor
## 0.3841799cor.test(s2_perf_overall$m, o_min_c_perf_s2)##
## Pearson's product-moment correlation
##
## data: s2_perf_overall$m and o_min_c_perf_s2
## t = -0.0051584, df = 136, p-value = 0.9959
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1675348 0.1666748
## sample estimates:
## cor
## -0.0004423248ComparePars(s2_perf_overall$m, o_min_c_perf_s2, use_identity_line = 0)
Invalid relationships to perf diffs
s1_o_inv <- s1_train %>% group_by(ID) %>% summarize(o_inv=mean(overt_invalid))
s1_c_inv <- s1_train %>% group_by(ID) %>% summarize(c_inv=mean(cognitive_invalid))
assert(all(s1_o_inv$ID==s1_c_inv$ID))
s1_inv_diff <- s1_c_inv$c_inv-s1_o_inv$o_inv
cor.test(o_min_c_perf_s1, s1_inv_diff)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and s1_inv_diff
## t = 0.029401, df = 123, p-value = 0.9766
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1730371 0.1781756
## sample estimates:
## cor
## 0.002650977ComparePars(o_min_c_perf_s1, s1_inv_diff, use_identity_line = 0)
s2_o_inv <- s2_train %>% group_by(ID) %>% summarize(o_inv=mean(prop_invalid_overt))
s2_c_inv <- s2_train %>% group_by(ID) %>% summarize(c_inv=mean(prop_invalid_cognitive))
assert(all(s2_o_inv$ID==s2_c_inv$ID))
s2_inv_diff <- s2_c_inv$c_inv-s2_o_inv$o_inv
cor.test(o_min_c_perf_s2, s2_inv_diff)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and s2_inv_diff
## t = 2.7091, df = 136, p-value = 0.007615
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.0614962 0.3790481
## sample estimates:
## cor
## 0.2262758ComparePars(s2_inv_diff, o_min_c_perf_s2, use_identity_line = 0, xchar="cog minus overt invalid", ychar="overt minus cog correct")
Percent of pts who actually showing better performance in cognitive
Training and test study 1
# Study 1
length(which(o_min_c_perf_s1 < 0))/length(o_min_c_perf_s1)## [1] 0.376length(which(o_min_c_perf_s1 == 0))/length(o_min_c_perf_s1) # 2.4% show no diff ## [1] 0.024length(which(o_min_c_testperf_s1 < 0))/length(o_min_c_testperf_s1)## [1] 0.328length(which(o_min_c_testperf_s1 == 0))/length(o_min_c_testperf_s1) ## [1] 0.104Training and test study 2
length(which(o_min_c_perf_s2 < 0))/length(o_min_c_perf_s2)## [1] 0.3913043length(which(o_min_c_perf_s2 == 0))/length(o_min_c_perf_s2) # < 1% show no diff## [1] 0.007246377length(which(o_min_c_testperf_s2 < 0))/length(o_min_c_testperf_s2)## [1] 0.3550725length(which(o_min_c_testperf_s2 == 0))/length(o_min_c_testperf_s2) ## [1] 0.1884058Modeling Results (Note: model numbers match those from
s.R but not those used in paper…)
… because the paper numbering follows the logic of how presented there, and because we culled poor/uninformative/auxiliary models as developing the streamlined final set
Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter
Numbers from before bug fix on the use of prior (3/23/23) dont use: 58986 | 52624 | 57510 | 74680
m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")
m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")
m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m1_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
# write.csv(m1_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")cat("\nBeta Median \n")##
## Beta Medianmedian(m1_study1_eb$beta)## [1] 2.330957cat("SD \n")## SDsd(m1_study1_eb$beta)## [1] 5.393476cat("\nRL LR \n")##
## RL LRmedian(m1_study1_eb$q_LR)## [1] 0.5916992cat("SD \n")## SDsd(m1_study1_eb$q_LR)## [1] 0.2752539cat("\nphi cog \n")##
## phi cogmedian(m1_study1_eb$cog_phi)## [1] 0.1593268cat("SD \n")## SDsd(m1_study1_eb$cog_phi)## [1] 0.2250586cat("\nphi overt \n")##
## phi overtmedian(m1_study1_eb$overt_phi)## [1] 0.1328623cat("SD \n")## SDsd(m1_study1_eb$overt_phi)## [1] 0.169769cat("\nepsilon \n")##
## epsilonmedian(m1_study1_eb$epsilon)## [1] 0.005554817cat("SD \n")## SDsd(m1_study1_eb$epsilon)## [1] 0.02699492cat("\nChoice LR \n")##
## Choice LRmedian(m1_study1_eb$choice_LR)## [1] 0.1198229cat("SD \n")## SDsd(m1_study1_eb$choice_LR)## [1] 0.1222895cat("\nES \n")##
## ESmedian(m1_study1_eb$explor_scalar)## [1] 0.3783385cat("SD \n")## SDsd(m1_study1_eb$explor_scalar)## [1] 0.2089571cat("\nBeta Median \n")##
## Beta Medianmedian(m1_study2_eb$beta)## [1] 2.692798cat("SD \n")## SDsd(m1_study2_eb$beta)## [1] 4.910113cat("\nRL LR \n")##
## RL LRmedian(m1_study2_eb$q_LR)## [1] 0.5543408cat("SD \n")## SDsd(m1_study2_eb$q_LR)## [1] 0.2806cat("\nphi cog \n")##
## phi cogmedian(m1_study2_eb$cog_phi)## [1] 0.2033696cat("SD \n")## SDsd(m1_study2_eb$cog_phi)## [1] 0.2858507cat("\nphi overt \n")##
## phi overtmedian(m1_study2_eb$overt_phi)## [1] 0.157541cat("SD \n")## SDsd(m1_study2_eb$overt_phi)## [1] 0.2223507cat("\nepsilon \n")##
## epsilonmedian(m1_study2_eb$epsilon)## [1] 7.844961e-07cat("SD \n")## SDsd(m1_study2_eb$epsilon)## [1] 0.02668559cat("\nChoice LR \n")##
## Choice LRmedian(m1_study2_eb$choice_LR)## [1] 0.09913484cat("SD \n")## SDsd(m1_study2_eb$choice_LR)## [1] 0.128331cat("\nES \n")##
## ESmedian(m1_study2_eb$explor_scalar)## [1] 0.2382697cat("SD \n")## SDsd(m1_study2_eb$explor_scalar)## [1] 0.1951073ComparePars(m1_study1_eb$cog_phi, m1_study1_eb$overt_phi,
"Phi", "Cog", "Overt")
ComparePars(m1_study2_eb$cog_phi, m1_study2_eb$overt_phi,
"Phi", "Cog", "Overt")
Relationships at both train and test in studies 1 and 2
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m1_study1_eb$ID)
o_min_c_phi_s1_m1 <- m1_study1_eb$overt_phi - m1_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m1
## t = -6.9629, df = 123, p-value = 1.765e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6469195 -0.3927853
## sample estimates:
## cor
## -0.5317171cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m1
## t = -6.0126, df = 123, p-value = 1.926e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6018399 -0.3284889
## sample estimates:
## cor
## -0.4766035assert(s2_perf[s2_perf$condition=="overt", "ID"]==m1_study2_eb$ID)
o_min_c_phi_s2_m1 <- m1_study2_eb$overt_phi - m1_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m1
## t = -4.7721, df = 136, p-value = 4.642e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5133414 -0.2259169
## sample estimates:
## cor
## -0.3787242cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m1
## t = -7.1401, df = 136, p-value = 5.093e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6339502 -0.3889985
## sample estimates:
## cor
## -0.5221611Descriptives and chi-square tests
Higher median decay in cog in both
# hist(sqrt(m1_study1_eb$cog_phi), breaks=100)
# hist(m1_study1_eb$cog_phi, breaks=100)
# hist(m1_study1_eb$overt_phi, breaks=100)
cat("\n Study 1 cog\n")##
## Study 1 cogmedian(m1_study1_eb$cog_phi)## [1] 0.1593268cat("\n Study 1 overt\n")##
## Study 1 overtmedian(m1_study1_eb$overt_phi)## [1] 0.1328623m1_s1_phi <- rbind(
data.frame("phi"=m1_study1_eb$overt_phi, "cond"="overt"),
data.frame("phi"=m1_study1_eb$cog_phi, "cond"="cognitive")
)
m1_s2_phi <- rbind(
data.frame("phi"=m1_study2_eb$overt_phi, "cond"="overt"),
data.frame("phi"=m1_study2_eb$cog_phi, "cond"="cognitive")
)
m1_s1_phi$Experiment <- "Experiment 1"
m1_s2_phi$Experiment <- "Experiment 2"
m1_phi_ID_comb <- rbind(m1_s1_phi, m1_s2_phi)
phi_medians <- rbind(
data.frame(m1_s1_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 1"),
data.frame(m1_s2_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 2")
)Difference higher for means but not as meaningful bc of skew, so plotting median
m1_s1_phi %>% group_by(cond) %>% summarize(m=mean(phi))## # A tibble: 2 × 2
## cond m
## <chr> <dbl>
## 1 cognitive 0.235
## 2 overt 0.176m1_s2_phi %>% group_by(cond) %>% summarize(m=mean(phi))## # A tibble: 2 × 2
## cond m
## <chr> <dbl>
## 1 cognitive 0.310
## 2 overt 0.209phi_plot <-ggplot(phi_medians, aes(x=cond, y=m, color=cond, group=cond)) +
geom_jitter(data = m1_phi_ID_comb,
aes(x=cond, y=phi, group=cond),
color="black", fill="gray57", size=2, pch=21, width=.15) +
geom_bar(stat="identity", fill="white", size=3, alpha=.2) +
facet_wrap(~ Experiment) + ga + ap + ft + lp + ylab(TeX('$\\phi') ) + xlab("") +
scale_color_manual(values=c("purple", "orange")) + ylim(-.03, 1.1) + tol
phi_plot
#ggsave("../paper/figs/fig_parts/phi_plot.png", phi_plot, height=5, width=10, dpi=700)Percent higher phi in studies 1 and 2
c_min_o_phi_s1_m1 <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phi
c_min_o_phi_s2_m1 <- m1_study2_eb$cog_phi - m1_study2_eb$overt_phitable((c_min_o_phi_s1_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s1_m1 > 0)*1))## 1
## 0.576table((c_min_o_phi_s2_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s2_m1 > 0)*1))## 1
## 0.6304348chisq.test(table((c_min_o_phi_s1_m1 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((c_min_o_phi_s1_m1 > 0) * 1)
## X-squared = 2.888, df = 1, p-value = 0.08924chisq.test(table((c_min_o_phi_s2_m1 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((c_min_o_phi_s2_m1 > 0) * 1)
## X-squared = 9.3913, df = 1, p-value = 0.00218m1_s1_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s1, "test_oc_diff"=o_min_c_testperf_s1, "phi_co_diff"=c_min_o_phi_s1_m1)
m1_s2_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s2, "test_oc_diff"=o_min_c_testperf_s2, "phi_co_diff"=c_min_o_phi_s2_m1)Training plots
a <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="black") +
geom_point(pch=21, fill="gray57", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("Proportion correct: \nOvert minus cognitive") +
xlab(TeX("")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
a## `geom_smooth()` using formula = 'y ~ x'
b <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="black") +
geom_point(pch=21, fill="gray57", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("") +
xlab(TeX("")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
b## `geom_smooth()` using formula = 'y ~ x'
Test plots
c <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
geom_point(pch=21, fill="gray79", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("Proportion correct: \nOvert minus cognitive") +
xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
d <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
geom_point(pch=21, fill="gray79", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("") +
xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)all_vs <- (a + b) / (c + d)all_vs## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#ggsave("../paper/figs/fig_parts/phi_vs_propdiff.png", all_vs, height=8, width=14, dpi=700)Correlation with delay effect
cor.test(m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi, ranef(m1_full_delay_cond)$ID$delay_cc1)##
## Pearson's product-moment correlation
##
## data: m1_study1_eb_v1$cog_phi + m1_study1_eb_v1$overt_phi and ranef(m1_full_delay_cond)$ID$delay_cc1
## t = -3.7767, df = 123, p-value = 0.0002462
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4712817 -0.1555501
## sample estimates:
## cor
## -0.3223523cor.test(m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi, ranef(m2_full_delay_cond)$ID$delay_cc1)##
## Pearson's product-moment correlation
##
## data: m1_study2_eb_v1$cog_phi + m1_study2_eb_v1$overt_phi and ranef(m2_full_delay_cond)$ID$delay_cc1
## t = -5.6966, df = 136, p-value = 7.246e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5646058 -0.2933200
## sample estimates:
## cor
## -0.4389115ComparePars((m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi), ranef(m1_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)
ComparePars((m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi), ranef(m2_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)
For resub, examining whether phi correlates with RT difference
assert(all(s1_tr_diffs_inds$ID==m1_study1_eb$ID))
assert(all(s1_te_diffs_inds$ID==m1_study1_eb$ID))
assert(all(s2_tr_diffs_inds$ID==m1_study2_eb$ID))
assert(all(s2_te_diffs_inds$ID==m1_study2_eb$ID))
ComparePars(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_tr_diffs_inds$m, use_identity_line = 0)
ComparePars(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_te_diffs_inds$m, use_identity_line = 0)
ComparePars(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_tr_diffs_inds$m, use_identity_line = 0)
ComparePars(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_te_diffs_inds$m, use_identity_line = 0)
cor.test(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_tr_diffs_inds$m)##
## Pearson's product-moment correlation
##
## data: m1_study1_eb$cog_phi - m1_study1_eb$overt_phi and s1_tr_diffs_inds$m
## t = 2.2364, df = 123, p-value = 0.02712
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.02286048 0.36075890
## sample estimates:
## cor
## 0.1976744cor.test(m1_study1_eb$cog_phi-m1_study1_eb$overt_phi, s1_te_diffs_inds$m)##
## Pearson's product-moment correlation
##
## data: m1_study1_eb$cog_phi - m1_study1_eb$overt_phi and s1_te_diffs_inds$m
## t = 3.2832, df = 123, p-value = 0.001336
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1139323 0.4376519
## sample estimates:
## cor
## 0.2838605cor.test(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_tr_diffs_inds$m)##
## Pearson's product-moment correlation
##
## data: m1_study2_eb$cog_phi - m1_study2_eb$overt_phi and s2_tr_diffs_inds$m
## t = 4.1166, df = 136, p-value = 6.634e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1755209 0.4736233
## sample estimates:
## cor
## 0.3328628cor.test(m1_study2_eb$cog_phi-m1_study2_eb$overt_phi, s2_te_diffs_inds$m)##
## Pearson's product-moment correlation
##
## data: m1_study2_eb$cog_phi - m1_study2_eb$overt_phi and s2_te_diffs_inds$m
## t = 2.6808, df = 136, p-value = 0.008252
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.05914673 0.37702705
## sample estimates:
## cor
## 0.2240372For resub, examining if phi differs based on condition order
unique_IDs <- unique(m1_study1_eb$ID)
for (s in 1:length(unique_IDs)) {
m1_study1_eb[m1_study1_eb$ID == unique_IDs[s], "first_cond"] <-
unique(s1_with_cond_order[s1_with_cond_order$deident_ID == unique_IDs[s], "first_cond"])
}unique_IDs <- unique(m1_study2_eb$ID)
for (s in 1:length(unique_IDs)) {
m1_study2_eb[m1_study2_eb$ID == unique_IDs[s], "first_cond"] <-
unique(s2_with_cond_order[s2_with_cond_order$de_ID == unique_IDs[s], "first_cond"])
}m1_study1_eb %>% group_by(first_cond) %>% summarize(om=median(overt_phi), cm=median(cog_phi))## # A tibble: 2 × 3
## first_cond om cm
## <chr> <dbl> <dbl>
## 1 cognTraining 0.118 0.180
## 2 overtTraining 0.148 0.133m1_study2_eb %>% group_by(first_cond) %>% summarize(om=median(overt_phi), cm=median(cog_phi))## # A tibble: 2 × 3
## first_cond om cm
## <chr> <dbl> <dbl>
## 1 cognitive 0.126 0.354
## 2 overt 0.173 0.141For resub, show delay distribution
hist(s1_train$delay, breaks=100)
hist(s2_train$delay, breaks=100)
Model validation, main text parts
6/5/23 — reran sims to make sure completely updated/bug fixed version of par results used
s1_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")
s1_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")a <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s1_train, s1_train_sim_m1_eb))
b <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s2_train, s2_train_sim_m1_eb))c <- PlotSimEmpTest(emp_df=s1_sit, sim_df=s1_test_sim_m1_eb)
d <- PlotSimEmpTest(emp_df=s2_sit, sim_df=s2_test_sim_m1_eb)a
b
sim_test_comb <- c + d #(temporarily flipped tol -> lp to get the legend for paper)
sim_test_comb
#ggsave("../paper/figs/fig_parts/train_sim_experiment1_plot.png", a, height=7, width=12, dpi=700)ggsave("../paper/figs/fig_parts/test_sim_experiment1_plot.png", c, height=5, width=11, dpi=700)ggsave("../paper/figs/fig_parts/test_sim_experiment2_plot.png", d, height=5, width=11, dpi=700)For supplemental because of space constraints
#ggsave("../paper/figs/fig_parts/train_sim_experiment2_plot.png", b, height=7, width=12, dpi=700)Model 3 — Same as m1 but not letting decay vary
Prior to vactoin — bug fixed files swapped in (see to_do notes)
m3_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior35433.csv")
m3_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior22212.csv")
m3_study1_eb <- rbind(m3_study1_eb_v1, m3_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m3_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior29818.csv")
m3_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior56742.csv")
m3_study2_eb <- rbind(m3_study2_eb_v1,m3_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m3_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")
# write.csv(m3_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")Worse in AIC
ComparePars(m3_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m3", "m1")
ComparePars(m3_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m3", "m1")
sum(m3_study1_eb$AIC-m1_study1_eb$AIC)## [1] 134.026sum(m3_study2_eb$AIC-m1_study2_eb$AIC)## [1] 147.2234Sanity check that the amount of model improvement when letting decay varies correlates with the diff in phi between conditions
assert(m1_study1_eb$ID==m3_study1_eb$ID)
assert(m1_study2_eb$ID==m3_study2_eb$ID)cor.test(m3_study1_eb$AIC-m1_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), method="spearman")##
## Spearman's rank correlation rho
##
## data: m3_study1_eb$AIC - m1_study1_eb$AIC and abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi)
## S = 107038, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6711582cor.test(m3_study2_eb$AIC-m1_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), method="spearman")##
## Spearman's rank correlation rho
##
## data: m3_study2_eb$AIC - m1_study2_eb$AIC and abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi)
## S = 163526, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6266436(Function prints Pearson’s r but more appropriate Spearman’s \(\rho\) reported one cell above)
ComparePars(m1_study1_eb$AIC-m3_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)
ComparePars(m1_study2_eb$AIC-m3_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)
Model 4 — Same as m3 but to 0 inits
m4_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits47402.csv")
m4_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits84391.csv")
m4_study1_eb <- rbind(m4_study1_eb_v1, m4_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m4_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits41957.csv")
m4_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits37660.csv")
m4_study2_eb <- rbind(m4_study2_eb_v1,
m4_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m4_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")
# write.csv(m4_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")Worse in AIC than pessimistic priors
ComparePars(m4_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m3", "m1")
ComparePars(m4_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m3", "m1")
sum(m4_study1_eb$AIC-m3_study1_eb$AIC)## [1] 120.8681sum(m4_study2_eb$AIC-m3_study2_eb$AIC)## [1] 446.3359Model 5 — Same as m1 but lets beta vary instead of decay
4/6 — bug fixed
m5_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta67330.csv")
m5_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta37484.csv")
m5_study1_eb <- rbind(m5_study1_eb_v1, m5_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m5_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta76602.csv")
m5_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta11129.csv")
m5_study2_eb <- rbind(m5_study2_eb_v1, m5_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Worse in AIC than diff decay
ComparePars(m5_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m5", "m1")
ComparePars(m5_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m5", "m1")
sum(m5_study1_eb$AIC-m1_study1_eb$AIC)## [1] 672.2797sum(m5_study2_eb$AIC-m1_study2_eb$AIC)## [1] 857.6129Model 6 — Same as m1 but lets learning rate vary instead of decay
4/6 — bug fix files swapped in
m6_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR12095.csv")
m6_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR13604.csv")
m6_study1_eb <- rbind(m6_study1_eb_v1, m6_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m6_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR69123.csv")
m6_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR43205.csv")
m6_study2_eb <- rbind(m6_study2_eb_v1, m6_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Worse in AIC than decay allowed to vary by condition
ComparePars(m6_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m6", "m1")
ComparePars(m6_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m6", "m1")
sum(m6_study1_eb$AIC-m1_study1_eb$AIC)## [1] 131.9173sum(m6_study2_eb$AIC-m1_study2_eb$AIC)## [1] 173.9311Higher overt learning rate in both
median(m6_study1_eb$cog_q_LR)## [1] 0.5048065median(m6_study1_eb$overt_q_LR)## [1] 0.556566median(m6_study2_eb$cog_q_LR)## [1] 0.4577678median(m6_study1_eb$overt_q_LR)## [1] 0.556566Learning rate correlated across conditions
ComparePars(m6_study1_eb$cog_q_LR, m6_study1_eb$overt_q_LR,
"q_LR", "Cog", "Overt")
ComparePars(m6_study2_eb$cog_q_LR, m6_study2_eb$overt_q_LR,
"q_LR", "Cog", "Overt")
Correlated with perf diffs but not as strongly as phi
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m6_study1_eb$ID)
o_min_c_q_LR_s1_m6 <- m6_study1_eb$overt_q_LR - m6_study1_eb$cog_q_LR
cor.test(o_min_c_perf_s1, o_min_c_q_LR_s1_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_q_LR_s1_m6
## t = 4.361, df = 123, p-value = 2.708e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2034049 0.5088491
## sample estimates:
## cor
## 0.3659413cor.test(o_min_c_testperf_s1, o_min_c_q_LR_s1_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_q_LR_s1_m6
## t = 2.6939, df = 123, p-value = 0.008049
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06303996 0.39525878
## sample estimates:
## cor
## 0.2360345assert(s2_perf[s2_perf$condition=="overt", "ID"]==m6_study2_eb$ID)
o_min_c_q_LR_s2_m6 <- m6_study2_eb$overt_q_LR - m6_study2_eb$cog_q_LR
cor.test(o_min_c_perf_s2, o_min_c_q_LR_s2_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_q_LR_s2_m6
## t = 3.83, df = 136, p-value = 0.0001949
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1528938 0.4553870
## sample estimates:
## cor
## 0.3120266cor.test(o_min_c_testperf_s2, o_min_c_q_LR_s2_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_q_LR_s2_m6
## t = 5.4029, df = 136, p-value = 2.849e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2724024 0.5489175
## sample estimates:
## cor
## 0.420372chisq.test(table((o_min_c_q_LR_s1_m6 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((o_min_c_q_LR_s1_m6 > 0) * 1)
## X-squared = 6.728, df = 1, p-value = 0.009491chisq.test(table((o_min_c_q_LR_s2_m6 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((o_min_c_q_LR_s2_m6 > 0) * 1)
## X-squared = 1.8551, df = 1, p-value = 0.1732Model 11 — Simple q-learner model
m11_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner24681.csv")
m11_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner77403.csv")
m11_study1_eb <- rbind(m11_study1_eb_v1,m11_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m11_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner18790.csv")
m11_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner26119.csv")
m11_study2_eb <- rbind(m11_study2_eb_v1, m11_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Much worse in AIC to not include decay
ComparePars(m11_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m11", "m3")
ComparePars(m11_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m11", "m3")
sum(m11_study1_eb$AIC-m3_study1_eb$AIC)## [1] 918.2861sum(m11_study2_eb$AIC-m3_study2_eb$AIC)## [1] 1225.945Model 12 — Simple Bayes model
m12_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner13415.csv")
m12_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner42048.csv")
m12_study1_eb <- rbind(m12_study1_eb_v1,m12_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m12_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner22977.csv")
m12_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner23536.csv")
m12_study2_eb <- rbind(m12_study2_eb_v1, m12_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))In turn much worse than the simple Q-learning model
ComparePars(m12_study1_eb$AIC, m11_study1_eb$AIC, "AIC", "m12", "m11")
ComparePars(m12_study2_eb$AIC, m11_study2_eb$AIC, "AIC", "m12", "m11")
sum(m12_study1_eb$AIC-m11_study1_eb$AIC)## [1] 622.5435sum(m12_study2_eb$AIC-m11_study2_eb$AIC)## [1] 1012.796Model 13 — M12 but allowing beta to vary
m13_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta60557.csv")
m13_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta85855.csv")
m13_study1_eb <- rbind(m13_study1_eb_v1, m13_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m13_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb <- rbind(m13_study2_eb_v1, m13_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Way worse than m1
ComparePars(m13_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m13", "m1")
ComparePars(m13_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m13", "m1")
sum(m13_study1_eb$AIC-m1_study1_eb$AIC)## [1] 2081.633sum(m13_study2_eb$AIC-m1_study2_eb$AIC)## [1] 2997.104Higher median beta in overt
median(m13_study1_eb$cog_beta)## [1] 1.790819median(m13_study1_eb$overt_beta)## [1] 1.876175median(m13_study2_eb$cog_beta)## [1] 1.893116median(m13_study2_eb$overt_beta)## [1] 2.099269Modest correlations between beta and perf
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m13_study1_eb$ID)
o_min_c_beta_s1_m13 <- m13_study1_eb$overt_beta - m13_study1_eb$cog_beta
cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_beta_s1_m13
## t = 3.5966, df = 123, p-value = 0.0004653
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1404791 0.4592086
## sample estimates:
## cor
## 0.3084768cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_beta_s1_m13
## t = 2.7331, df = 123, p-value = 0.007199
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06646281 0.39815509
## sample estimates:
## cor
## 0.2392776assert(s2_perf[s2_perf$condition=="overt", "ID"]==m13_study2_eb$ID)
o_min_c_beta_s2_m13 <- m13_study2_eb$overt_beta - m13_study2_eb$cog_beta
cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_beta_s2_m13
## t = 2.8955, df = 136, p-value = 0.004412
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.07696563 0.39227998
## sample estimates:
## cor
## 0.2409713cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_beta_s2_m13
## t = 1.9265, df = 136, p-value = 0.05613
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.004233115 0.321339865
## sample estimates:
## cor
## 0.162987Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond
Alt refers to this being the simpler version with uncertainty factored in at the action (arm) level rather than a weighted uncertainty computed — this and the more complex version were v similar in terms of fit
m14_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt39650.csv")
m14_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt71571.csv")
m14_study1_eb <- rbind(m14_study1_eb_v1, m14_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m14_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt38226.csv")
m14_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt28286.csv")
m14_study2_eb <- rbind(m14_study2_eb_v1, m14_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Original version with the more complicated weighting scheme
m14_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes25885.csv")
m14_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes65430.csv")
m14_study1_eb_old <- rbind(m14_study1_eb_v1_old, m14_study1_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))
m14_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes30853.csv")
m14_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes58630.csv")
m14_study2_eb_old <- rbind(m14_study2_eb_v1_old, m14_study2_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))As expected, these are very similar, with the new version a tad better
ComparePars(m14_study1_eb$AIC, m14_study1_eb_old$AIC, "AIC", "m14", "m14 old")
sum(m14_study1_eb$AIC-m14_study1_eb_old$AIC)## [1] -18.67696ComparePars(m14_study2_eb$AIC, m14_study2_eb_old$AIC, "AIC", "m14", "m14 old")
sum(m14_study2_eb$AIC-m14_study2_eb_old$AIC)## [1] -43.19563Worse AIC (and nll) compared to m1 (not a nested model bc doesn’t have diff decay)
ComparePars(m14_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m14", "m1")
ComparePars(m14_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m14", "m1")
sum(m14_study1_eb$AIC-m1_study1_eb$AIC)## [1] 448.2412sum(m14_study2_eb$AIC-m1_study2_eb$AIC)## [1] 486.6508sum(m14_study1_eb$nll-m1_study1_eb$nll)## [1] 99.12058sum(m14_study2_eb$nll-m1_study2_eb$nll)## [1] 105.3254Higher median rho in overt
median(m14_study1_eb$rho_cog)## [1] 0.1651074median(m14_study1_eb$rho_overt)## [1] 0.2611568median(m14_study2_eb$rho_cog)## [1] 0.1066345median(m14_study2_eb$rho_overt)## [1] 0.2367347Correlations w perf but not as strongly as phi
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m14_study1_eb$ID)
o_min_c_beta_s1_m14 <- m14_study1_eb$rho_overt - m14_study1_eb$rho_cog
cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_beta_s1_m14
## t = 4.2195, df = 123, p-value = 4.71e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1919767 0.4999821
## sample estimates:
## cor
## 0.3555962cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_beta_s1_m14
## t = 3.3007, df = 123, p-value = 0.001262
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1154251 0.4388741
## sample estimates:
## cor
## 0.2852507assert(s2_perf[s2_perf$condition=="overt", "ID"]==m14_study2_eb$ID)
o_min_c_beta_s2_m14 <- m14_study2_eb$rho_overt - m14_study2_eb$rho_cog
cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_beta_s2_m14
## t = 5.1757, df = 136, p-value = 7.977e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2558970 0.5363993
## sample estimates:
## cor
## 0.4056554cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_beta_s2_m14
## t = 3.4613, df = 136, p-value = 0.0007187
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1232954 0.4311429
## sample estimates:
## cor
## 0.2845379Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond
Alternate version with the simpler action-wise weighting (see m14 description)
m15_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt33317.csv")
m15_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt35801.csv")
m15_study1_eb <- rbind(m15_study1_eb_v1, m15_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m15_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt64665.csv")
m15_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt54458.csv")
m15_study2_eb <- rbind(m15_study2_eb_v1, m15_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Original version with the more complicated weighting scheme
bug numbers: 75340 81994 63931 40889
m15_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes50404.csv")
m15_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes56900.csv")
m15_study1_eb_old <- rbind(m15_study1_eb_v1_old, m15_study1_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))
m15_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes68726.csv")
m15_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes30103.csv")
m15_study2_eb_old <- rbind(m15_study2_eb_v1_old, m15_study2_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))As expected, again very similar — so using new/simpler one going forward
ComparePars(m15_study1_eb$AIC, m15_study1_eb_old$AIC, "AIC", "m15", "m15 old")
sum(m15_study1_eb$AIC-m15_study1_eb_old$AIC)## [1] 11.82675ComparePars(m15_study2_eb$AIC, m15_study2_eb_old$AIC, "AIC", "m15", "m15 old")
sum(m15_study2_eb$AIC-m15_study2_eb_old$AIC)## [1] -53.60774Worse AIC compared to m1
ComparePars(m15_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m15", "m1")
ComparePars(m15_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m15", "m1")
sum(m15_study1_eb$AIC-m1_study1_eb$AIC)## [1] 169.3991sum(m15_study2_eb$AIC-m1_study2_eb$AIC)## [1] 127.4026# sum(m15_study1_eb$nll-m1_study1_eb$nll)
# sum(m15_study2_eb$nll-m1_study2_eb$nll)Higher median decay in cog in both
median(m15_study1_eb$cog_phi)## [1] 0.2204569median(m15_study1_eb$overt_phi)## [1] 0.1742236median(m15_study2_eb$cog_phi)## [1] 0.3308409median(m15_study1_eb$overt_phi)## [1] 0.1742236ComparePars(m15_study1_eb$cog_phi, m15_study1_eb$overt_phi,
"Phi", "Cog", "Overt")
ComparePars(m15_study2_eb$cog_phi, m15_study2_eb$overt_phi,
"Phi", "Cog", "Overt")
Similar perf correlations as m1
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m15_study1_eb$ID)
o_min_c_phi_s1_m15 <- m15_study1_eb$overt_phi - m15_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m15
## t = -7.1705, df = 123, p-value = 6.097e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6560039 -0.4060503
## sample estimates:
## cor
## -0.542943cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m15
## t = -6.1734, df = 123, p-value = 8.902e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6098839 -0.3397791
## sample estimates:
## cor
## -0.4863663assert(s2_perf[s2_perf$condition=="overt", "ID"]==m15_study2_eb$ID)
o_min_c_phi_s2_m15 <- m15_study2_eb$overt_phi - m15_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m15
## t = -3.8478, df = 136, p-value = 0.0001825
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4565342 -0.1543079
## sample estimates:
## cor
## -0.3133333cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m15
## t = -6.5236, df = 136, p-value = 1.247e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6058791 -0.3496199
## sample estimates:
## cor
## -0.4882024Model 27 — Fixing both ES and epsilon at the group level
m27_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed46645.csv")
m27_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed33416.csv")
m27_study1_eb <- rbind(m27_study1_eb_v1, m27_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m27_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed20328.csv")
m27_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed44601.csv")
m27_study2_eb <- rbind(m27_study2_eb_v1, m27_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m27_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")
# write.csv(m27_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")Does worsen NLL substantially
ComparePars(m27_study1_eb$nll, m1_study1_eb$nll, "nll", "m27", "m1")
ComparePars(m27_study2_eb$nll, m1_study2_eb$nll, "nll", "m27", "m1")
sum(m27_study1_eb$nll-m1_study1_eb$nll)## [1] 259.0197sum(m27_study2_eb$nll-m1_study2_eb$nll)## [1] 367.588Still higher median decay in cog in both
median(m27_study1_eb$cog_phi)## [1] 0.1637271median(m27_study1_eb$overt_phi)## [1] 0.1337206median(m27_study2_eb$cog_phi)## [1] 0.2304866median(m27_study1_eb$overt_phi)## [1] 0.1337206Strong relationships at both train and test in studies 1 and 2
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m27_study1_eb$ID)
o_min_c_phi_s1_m27 <- m27_study1_eb$overt_phi - m27_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m27
## t = -6.9726, df = 123, p-value = 1.68e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6473518 -0.3934142
## sample estimates:
## cor
## -0.5322504cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m27
## t = -6.4826, df = 123, p-value = 1.966e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6248649 -0.3610156
## sample estimates:
## cor
## -0.504631assert(s2_perf[s2_perf$condition=="overt", "ID"]==m27_study2_eb$ID)
o_min_c_phi_s2_m27 <- m27_study2_eb$overt_phi - m27_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m27
## t = -4.8922, df = 136, p-value = 2.775e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5203133 -0.2349247
## sample estimates:
## cor
## -0.3868433cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m27
## t = -7.114, df = 136, p-value = 5.847e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6328069 -0.3873771
## sample estimates:
## cor
## -0.52077159 and 62% of pts show a higher phi in overt
table((o_min_c_phi_s1_m27 > 0)*1)##
## 0 1
## 73 52table((o_min_c_phi_s2_m27 > 0)*1)##
## 0 1
## 86 52table((o_min_c_phi_s1_m27 > 0))[1]/sum(table((o_min_c_phi_s1_m27 > 0)))## FALSE
## 0.584table((o_min_c_phi_s2_m27 > 0))[1]/sum(table((o_min_c_phi_s2_m27 > 0)))## FALSE
## 0.6231884chisq.test(table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0)))##
## Chi-squared test for given probabilities
##
## data: table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0))
## X-squared = 3.528, df = 1, p-value = 0.06034chisq.test(table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0)))##
## Chi-squared test for given probabilities
##
## data: table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0))
## X-squared = 8.3768, df = 1, p-value = 0.0038Model 29 — Model 1 but with no choice kernel
m29_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK21030.csv")
m29_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK89190.csv")
# Replace when back
#m29_study1_eb <- m29_study1_eb_v1
m29_study1_eb <- rbind(m29_study1_eb_v1, m29_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m29_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK64103.csv")
m29_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK74167.csv")
# Delete when final one is back
m29_study2_eb <- m29_study2_eb_v1
m29_study2_eb <- rbind(m29_study2_eb_v1, m29_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m29_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")
# write.csv(m29_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")Dramatically worse fit
ComparePars(m29_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m29", "m1")
ComparePars(m29_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m29", "m1")
sum(m29_study1_eb$AIC-m1_study1_eb$AIC)## [1] 1815.325sum(m29_study2_eb$AIC-m1_study2_eb$AIC)## [1] 2777.038Model Comparison
AIC_diff_df <-
rbind(
data.frame("diff"=m12_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes"),
data.frame("diff"=m12_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes"),
data.frame("diff"=m4_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decayto0"),
data.frame("diff"=m4_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decayto0"),
data.frame("diff"=m3_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay"),
data.frame("diff"=m3_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay"),
data.frame("diff"=m5_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_beta"),
data.frame("diff"=m5_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_beta"),
data.frame("diff"=m13_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes_diff_beta"),
data.frame("diff"=m13_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes_diff_beta"),
data.frame("diff"=m6_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_lr"),
data.frame("diff"=m6_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_lr"),
data.frame("diff"=m14_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_bayes"),
data.frame("diff"=m14_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_bayes"),
data.frame("diff"=m15_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_diff_decay_bayes"),
data.frame("diff"=m15_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_diff_decay_bayes"),
data.frame("diff"=m1_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_phi"),
data.frame("diff"=m1_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_phi")
)
AIC_diff_df$study <- as.factor(AIC_diff_df$study)
AIC_diff_summ <- AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff))## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.sum(m1_study1_eb$AIC-m3_study1_eb$AIC)## [1] -134.026sum(m1_study2_eb$AIC-m3_study2_eb$AIC)## [1] -147.2234Confirm min in both studies is diff phi
AIC_diff_summ[AIC_diff_summ$study == 1, ] %>% arrange(m)## # A tibble: 9 × 3
## # Groups: model [9]
## model study m
## <chr> <fct> <dbl>
## 1 Q_decay_diff_phi 1 -8.42
## 2 Q_decay_diff_lr 1 -7.36
## 3 Q_decay 1 -7.35
## 4 Q_diff_decay_bayes 1 -7.06
## 5 Q_decayto0 1 -6.38
## 6 Q_decay_diff_bayes 1 -4.83
## 7 Q_decay_diff_beta 1 -3.04
## 8 Bayes 1 4.98
## 9 Bayes_diff_beta 1 8.23AIC_diff_summ[AIC_diff_summ$study == 2, ] %>% arrange(m)## # A tibble: 9 × 3
## # Groups: model [9]
## model study m
## <chr> <fct> <dbl>
## 1 Q_decay_diff_phi 2 -9.95
## 2 Q_diff_decay_bayes 2 -9.03
## 3 Q_decay 2 -8.88
## 4 Q_decay_diff_lr 2 -8.69
## 5 Q_decay_diff_bayes 2 -6.42
## 6 Q_decayto0 2 -5.65
## 7 Q_decay_diff_beta 2 -3.74
## 8 Bayes 2 7.34
## 9 Bayes_diff_beta 2 11.8levels <-
c(
"Bayes_diff_beta",
"Bayes",
"Q_decayto0",
"Q_decay",
"Q_decay_diff_beta",
"Q_decay_diff_lr",
"Q_decay_diff_bayes",
"Q_diff_decay_bayes",
"Q_decay_diff_phi"
)
assert(length(levels)==9)
AIC_diff_summ$model <- factor(AIC_diff_summ$model, levels = levels)
AIC_diff_df$model <- factor(AIC_diff_df$model, levels = levels)mp1 <- ggplot(AIC_diff_summ %>% filter(study == 1), aes(x=model, y=m, group=model, fill=model)) +
geom_hline(yintercept = 0, size=4) +
geom_hline(yintercept = seq(-50, 60, 10), color="gray57", alpha=.8, size=.65) +
geom_jitter(data=AIC_diff_df %>% filter(study == 1), aes(x=model, y=diff, group=model),
fill="gray57", pch=21, width=.3, size=4, alpha=.08) +
geom_bar(stat="identity", color="black", alpha=.8, size=2) +
ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
tol + ga + ap + tp +
ggtitle("Study 1") +
#facet_wrap(~ study) +
xlab("") +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
theme(axis.title.y=element_text(size=22)) + ylim(-55, 65) + coord_flip()mp2 <- ggplot(AIC_diff_summ %>% filter(study == 2), aes(x=model, y=m, group=model, fill=model)) +
geom_hline(yintercept = 0, size=4) +
geom_hline(yintercept = seq(-50, 60, 10), color="gray57", alpha=.8, size=.65) +
geom_jitter(data=AIC_diff_df %>% filter(study == 2), aes(x=model, y=diff, group=model),
fill="gray57", pch=21, width=.3, size=4, alpha=.08) +
geom_bar(stat="identity", color="black", alpha=.8, size=2) +
ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
tol + ga + ap + tp +
ggtitle("Study 2") +
xlab("") +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
theme(axis.title.y=element_text(size=22)) + ylim(-55, 65) + coord_flip()combined_aic <- mp1 + mp2combined_aic
#ggsave("../paper/figs/fig_parts/aic_plot.png", combined_aic, height=5, width=10, dpi=700)AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff)) %>% arrange(m)## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.## # A tibble: 18 × 3
## # Groups: model [9]
## model study m
## <fct> <fct> <dbl>
## 1 Q_decay_diff_phi 2 -9.95
## 2 Q_diff_decay_bayes 2 -9.03
## 3 Q_decay 2 -8.88
## 4 Q_decay_diff_lr 2 -8.69
## 5 Q_decay_diff_phi 1 -8.42
## 6 Q_decay_diff_lr 1 -7.36
## 7 Q_decay 1 -7.35
## 8 Q_diff_decay_bayes 1 -7.06
## 9 Q_decay_diff_bayes 2 -6.42
## 10 Q_decayto0 1 -6.38
## 11 Q_decayto0 2 -5.65
## 12 Q_decay_diff_bayes 1 -4.83
## 13 Q_decay_diff_beta 2 -3.74
## 14 Q_decay_diff_beta 1 -3.04
## 15 Bayes 1 4.98
## 16 Bayes 2 7.34
## 17 Bayes_diff_beta 1 8.23
## 18 Bayes_diff_beta 2 11.8Get AICs for supplemental table
cat(round(sum(m1_study1_eb$AIC), 0), "\n")## 42125cat(round(sum(m1_study2_eb$AIC), 0))## 42798cat(round(sum(m3_study1_eb$AIC), 0), "\n")## 42259cat(round(sum(m3_study2_eb$AIC), 0))## 42945cat(round(sum(m4_study1_eb$AIC), 0), "\n")## 42380cat(round(sum(m4_study2_eb$AIC), 0))## 43391cat(round(sum(m5_study1_eb$AIC), 0), "\n")## 42797cat(round(sum(m5_study2_eb$AIC), 0))## 43655cat(round(sum(m6_study1_eb$AIC), 0), "\n")## 42257cat(round(sum(m6_study2_eb$AIC), 0))## 42972cat(round(sum(m13_study1_eb$AIC), 0), "\n")## 44206cat(round(sum(m13_study2_eb$AIC), 0))## 45795cat(round(sum(m14_study1_eb$AIC), 0), "\n")## 42573cat(round(sum(m14_study2_eb$AIC), 0))## 43284cat(round(sum(m15_study1_eb$AIC), 0), "\n")## 42294cat(round(sum(m15_study2_eb$AIC), 0))## 42925cat(round(sum(m11_study1_eb$AIC), 0), "\n")## 43177cat(round(sum(m11_study2_eb$AIC), 0))## 44171cat(round(sum(m12_study1_eb$AIC), 0), "\n")## 43800cat(round(sum(m12_study2_eb$AIC), 0))## 45184